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Abstract: 

Gravitational light deflection can distort the images of distant sources by its tidal effects. 
The population of faint blue galaxies is at sufficiently high redshift so that their images 
are distorted near foreground clusters, with giant luminous arcs being the most spectacular 
evidence for this effect. Much weaker distortions, however, can observationally be detected 
by a statistical analysis of the numerous faint galaxy images, as first demonstrated by 
Tyson, Valdes & Wenk. This distortion effect can be used as a quantitative tool for the 
reconstruction of the surface mass density of galaxy clusters with appropriate redshifts, as 
was demonstrated by Kaiser & Squires. They have derived an explicit equation for this 
surface mass density in terms of its tidal field. 

The reconstruction formula by Kaiser & Squires must be modified because of two effects: 
in its original form it applies only to weak lenses, and hence must be generalized to account 
for stronger lensing effects. Second, due to the nature of the inversion formula, it produces 
boundary artefacts (or biases) if applied to real data which are confined to a finite field on 
the sky. We discuss several possibilities to obtain inversion formulae which are exact for ideal 
data on a finite data field (CCD). We demonstrate that there exists an infinite number of such 
finite-field inversion formulae, which differ in their sensitivity to observational effects (such as 
noise, intrinsic ellipticities of the sources, etc.). We show that, using two simple conditions, 
one can uniquely specify a finite-field formula which in a well-defined sense minimizes the 
sensitivity to observational effects. We then use synthetic data to compare the quality of 
our new reconstruction method with that of previous finite-field inversion techniques and of 
the nonlinear generalization of the Kaiser & Squires method. This analysis demonstrates 
that our new inversion method is superior to the other previously considered finite-field 
inversions, and only slightly more noisy than the Kaiser & Squires inversion, but in contrast 
to the latter, it lacks the boundary artefacts. We shall discuss that the lack of boundary 
artefacts and the slightly increased noise have the same origin, and that every finite-field 
reconstruction must be more noisy than that obtained from the Kaiser & Squires method. 
Furthermore, the rms deviations of the reconstructed density field from the input surface 
mass density are fairly homogeneously distributed over the field. We therefore conclude that 
the new method developed here is the best inversion formula yet found. 
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1 Introduction 

Weak distortions of images of faint, high-redshift galaxies due to the gravitational light 
deflection caused by a foreground cluster can be used as a probe for the mass distribution 
of the cluster (e.g., Webster 1985). This effect manifests itself most visibly in the giant 
luminous arcs, where the image distortion is very strong (see the recent review by Fort & 
Mellier 1994). It is, of course, more likely that a galaxy image is less dramatically distorted, 
and the first of those more weakly distorted images (called arclets) were found by Fort et al. 
(1988). With one giant luminous arc in A370 and several arclets, the mass distribution of 
the cluster could be much better constrained than from single arcs alone (e.g., Grossman & 
Narayan 1989, Kneib et al 1993). Even weaker distortions, not obviously seen in individual 
galaxy images but clearly detectable statistically were verified by Tyson, Valdes & Wenk 
(1990). Such observations can probe the cluster mass distribution to much larger angular 
scales from the cluster center than accessible with giant luminous arcs (Bonnet, Mellier & 
Fort 1994). 

The distortion of images can be calculated, loosely speaking, as a convolution of the 
surface mass density with a (known) weighting funtion. In their pioneering paper, Kaiser & 
Squires (1993, hereafter KS) have shown that this relation can be inverted, i.e., the surface 
mass density distribution can be obtained as a convolution of the distortion field with a 
(known) kernel (earlier qualitative work on the determination of cluster mass distributions 
from weak distortions include Kochanek 1990; Miralda-Escude 1991). In other words, if the 
distortion field can be measured with sufficient accuracy (i.e., if the observations extend 
to sufficiently faint limits, so that the number density of faint galaxies is high enough), 
a density map of the corresponding cluster can be constructed. In fact, this method has 
already been applied to a number of clusters (Fahlman et al. 1994; Smail et al. 1995; Kaiser 
et al. 1994, henceforth KSFWB), and in particular for the cluster MS1224 has led to a 
mass estimate which is significantly larger than estimates from a dynamical study of this 
cluster (Carlberg, Yee & Ellingson 1995). These first applications have amply demonstrated 
the potential power of this new method for investigating clusters, which is independent of 
assumptions about the physical or dynamical state of the cluster galaxies or the intracluster 
gas. 

Nevertheless, the application of the KS method is far from being trivial. The difficulties 
with it can roughly be split into two parts: first, and certainly most serious, is the problem 
of extracting the distortion field from observations. It is essential for the method that the 
observations are taken at very good seeing, since the typical size of the galaxies - at the 
faint levels for which their number density is sufficiently large - is comparable to the seeing 
disk, so that the distortion signal is diluted. Other effects, like telescope tracking errors, 
lead to a non-trivial point spread function whose anisotropy generates an artifical distortion. 
These problems can be overcome, as discussed in detail by Bonnet & Mellier (1995), Kaiser, 
Squires & Broadhurst (1995) and Mould et al. (1994), and are not the subject of this paper. 

The second set of problems is more of a theoretical nature. The KS inversion method 
relates the 'shear' of the mass distribution - which depends linearly on the surface mass den- 
sity - to the mass distribution. However, the shear is not an observable. As was pointed out 
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in Schneider & Seitz (1995, hereafter Paper I), the observable is a combination of the shear 
and the surface mass density; only in the weak-lensing regime can the shear be obtained from 
observable quantities. This fact complicates the inversion process due to its nonlinearity. In 
Seitz & Schneider (1995, hereafter Paper II) it was shown that a general two-dimensional 
nonlinear reconstruction is possible, so that the KS method can be extended towards the 
inner parts of the lensing cluster. This method, though slightly more complicated than the 
original KS method, works as well as the KS inversion in those regimes where both can 
be applied. A second problem with the KS method is the fact that the inversion formula 
in principle requires data on the entire sky, whereas the observed field around a cluster is 
always bounded. Hence, the shear field has to be extrapolated outside the data field. The 
simplest such 'extrapolation' is setting the shear to zero outside the data field. This, how- 
ever, leads to serious 'boundary effects' in the reconstructed density field. Alternatively, one 
can fit a model function to the shear field in the outer parts of the data field and use this 
model function for the extrapolation, as has been done in Paper II. Using a sample of cluster 
models drawn from a high-resolution CDM simulation, it was shown by Bartelmann (1995) 
that this extrapolation yields in fact fairly accurate estimates of the cluster masses, and the 
reconstructed density fields in Paper II show that the boundary artefacts obtained from the 
original KS inversion can be reduced significantly by this extrapolation. Nevertheless, both 
of these two extrapolations mean 'inventing data' and should thus be avoided. 

It was shown by Schneider (1995, hereafter Paper III), using a relation between the 
gradient of the surface mass density and a combination of derivatives of the shear which 
was obtained by Kaiser (1995), that a cluster inversion formula can be obtained which 
explicitly makes use only of data on a finite field. This formula was obtained by a carefully 
chosen average of line integrals of the gradient of the surface mass density over the data 
field. Subsequently, KSFWB have generalized the ideas of Paper III to obtain different sets 
of such finite-field inversion equations, one of which was explored in detail by Bartelmann 
(1995). We shall summarize these finite-field methods in Sect. 3 below, where it will become 
clear that one has a huge freedom for constructing such inversion formulae. It thus remains 
questionable which of these is the 'best' one - all of them yield exact results for perfect data, 
but they behave differently with respect to noise, both from observations and from the fact 
that the faint galaxies have an intrinsic ellipticity distribution. 

In this paper we shall reconsider the problem of finding a finite-field inversion equation. 
The notation and the nature of the problem are introduced in Sect. 2. In Sect. 3 we derive 
families of finite-field inversions based on 'averaging line integrals' of the gradient of the 
surface mass density, of which the ones of Paper III and of Bartelmann (1995) are special 
cases; this section summarizes the previous attempts for constructing finite-field inversion 
methods and generalizes them slightly. In Sect. 4 we derive a finite-field inversion formula 
which is uniquely specified by two well-motivated requirements. In Sect. 5 we compare the 
various inversion techniques (KS, the one of Paper III, that from Bartelmann 1995, and our 
new one) using simulated distortion fields. It turns out that our new inversion method is 
superior to the other finite-field methods, and that it has (at worst) only marginally worse 
noise properties than the KS method, but the advantage that it is free from (systematic) 
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boundary artefacts. The relative noise properties of our new method and KS depend strongly 
on the underlying mass distribution; if it extends to near the boundary of the data field, 
the noise properties of our inversion are better than that of KS. Also, we show that the 
deviations of our reconstructed mass profile from the original one are statistically homo- 
geneously distributed over the data field, in contrast to the other reconstruction formulae. 
We thus conclude that our new inversion technique is the best found up to now, since the 
slightly increased noise relative to KS in some situations is well compensated by the lack of 
systematic boundary artefacts. 

Throughout this paper we consider the possibility that the lens is not weak, so the linear 
approximation of the shear is not used. Hence, if we apply the KS method, we always do 
so by the iterative procedure developed in Paper II. However, we assume that the lenses are 
non-critical, so that the local degeneracy found in Paper I does not occur. Furthermore, as in 
the previous papers, we will assume that all sources are at the same redshift; this assumption 
simplifies the inversion process considerably, but is not essential for the application of our 
inversion technique. We shall deal with a redshift distribution of the sources in a later 
paper. If we were to consider only weak lenses, and thus use the linear approximation, then 
a redshift distribution of the sources causes no additional technical problems. 



2 Inversion methods 

We use the same notation as in Papers I— III. Briefly, we define the deflection potential 

ip{9) = - f dV k(0') In 1 0-0' I (2.1) 

in terms of the dimensionless surface mass density k{0). The linearized lens mapping /3 = 
— j \/i/j(0), which describes the distortion of small sources, is d/3 = A(0) d0, where (3 denotes 
the angular position on the source sphere, and 

A(9)=( 1 - K + ^ 1 + 72 ) (2.2) 

is the Jacobian matrix of the lens mapping, which is related to the deflection potential ip 
through 

K{9) = \v 2 m , 

7l(0) = ^,22-V>ll) , (2 ' 3) 
72 (0) = _^ jl2 , 

where indices separated by a comma denote partial derivatives with respect to 0;. For details 
concerning these lensing relations, see Schneider, Ehlers & Falco (1992). Combining (2.1 & 
2.3), and defining the complex shear 7(0) = 71 (9) + i7 2 (0), one obtains 
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7 (0) = - f d 2 0' V{0 - 0') k(0') , (2.4) 
7f Jr 2 



where 

\e 



V (°) = ~ 2 J 4 = (*) + ^2 W (2.5) 



is a complex kernel. It was shown in KS that (2.5) can be inverted to yield 

K (0) = - f d 2 6' ne\V*(0-0')^(0')} . (2.6a) 

T Jr 2 

Using partial integration and the assumption that the shear vanishes at infinity yields a 
different form of the inversion integral, 

^=^H^,,)(^[Jj;^) , (2.6.) 

with 

H KS (fl',fl) = — = V > (-— \n\0-0'\) ■ (2.6c) 

1 ' ^\0-0'\ 2 6 v 2tt 1 i; 1 ; 

hence, the surface mass density k at a position can be obtained by convolving the deflection 
angle of a point mass with negative mass —1/2 positioned at with the derivative of the 
shear field 7. 

As mentioned in the introduction, there are two fundamental problems associated with 
(2.6a). First, the shear 7(0) cannot be obtained from image distortions, but what can be 
determined is the (complex) distortion 

27(1 ; (2-7) 



(1 -k) 2 + | 7 | 

see Paper I and Miralda-Escude 1991 for details. If we assume (as we shall do in the rest of 
the paper) that the cluster is non-critical (i.e., detA(0) > for all 0), then the quantity 
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\5\ 2 



(2.8) 



(1-k) \^ \5\ 2 
is also an observable. Combining (2.8) with (2.6a), we obtain 

K (0) = - ( d 2 6' [1-K(0')]ne [V*(0-0')g{0')\ , (2.9) 

i.e., an integral equation for k(0), which can easily be solved iteratively, as shown in Paper II. 
A second problem associated with (2.6) and (2.9) is the range of integration: since data are 
available only on a finite field, given by the size of the CCD, these equations can only be 
applied to data if a 'guess' is made for g outside the data field U (where U is that region 
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within which we can determine g from observations) . Typically, g and thus 7 is set to zero 
outside the data field 1 , so that 

k ks,1(0) . = I f d 2 e ' ne [£>*(0-0') 7 (0')] . (2.10a) . 

n Ju 

This so called 'KS'-estimate suffers from boundary artefacts, whose amplitude depends on 
form and size of the data field and on the mass distribution. Alternatively, one can try to ex- 
trapolate g outside the field, which is reasonable since the shear cannot decrease faster than 
the shear of a point mass. Such an extrapolation was successfully applied to a numerical 
model cluster in Paper II, and in a more systematic study Bartelmann (1995) has demon- 
strated that the extrapolation method yields quite accurate results for the mass within 
circular apertures. One can obtain a different estimate for n from (2.6b) by setting the 
derivatives of 7 to zero outside the data field 

K ^(«):=£dVH-(«',«)(:^|^^^) . (2.105) 

This estimate is expected to differ from that in (2.10a) in the systematic boundary artefacts 
as well as in the noise properties. 

In order to determine mass profiles, one should aim for an inversion of (2.4) which 
makes use only of the observed data on U and which yields an exact reconstruction for 
perfect data. As was shown by KS, the surface mass density can be determined only up to 
an additive constant, if the shears is assumed to be measurable. 2 Choosing this constant to 
be the (unknown) mean surface mass density R over the data field, one may want to search 
for a kernel V(0' , 0) such that 

K (0)-R= - f d 2 6' ne\V*(0', 0)^(0')] . (2.11) 

7T Ju 

If we use (2.4) for 7 in the preceding equation, we see that the kernel V has to satisfy the 
relation 

Q(0, tf) : = J_ / d 2 6' Tie \f>(0', 0)V(0' - 0)1 = 5(0 - -0) - ,0eW,i5el 2 (2.12) 

TT Ju A 

where &(■&) = 1 if ■d e U, and zero otherwise, and A is the area enclosed by U. However, 
we note that (2.11) is not the only possible form of an inversion equation; a more general 
form would include an additional integral over the boundary of U. 

In the next two sections, we present two different strategies to obtain inversion equations 
which are exact on a finite data field. By that we mean that if 7 obeys (2.4), or, alternatively, 



■ 1 2 11 

1 Note that |-y| decreases at most as 1/ , and for an isothermal distribution, 7 oc 1/ . 

It is important to note that the shear 7 is not an observable in general; only for the case of weak lensing, 
the observable g - see Eq. (2.8) - can be linearized, and in this linear approximation, 7 becomes the 
observable. 
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if g obeys (2.8), then the surface mass density n can be obtained up to a constant. However, 
since the value of g has to be determined by averaging over observed galaxy images, the 
relation between (the observed) g and (the true) n will deviate from (2.8), due to noise 
coming from the discreteness of galaxy images, the intrinsic ellipticity of sources etc. In 
Sect. 3 we will summarize and generalize the work previously done on finite-field inversion 
formulae which are all exact for 'perfect data', but differ from their behaviour with respect 
to observational 'noise'. In Sect. 4, we will then obtain a new inversion method which is 
motivated by the fact that a certain 'noise component' can be identified as such and can 
be 'filtered out' in the reconstruction. We shall then, in Sect. 5, demonstrate that our new 
inversion formulae is superiour to the other finite-field reconstructions published (Paper III, 
KSFWB, Bartelmann 1995). 

3 Finite field kernels: Averaging over line integrals 

The first inversion formula on a finite field was derived in Paper III; its derivation started 
from the equation 

Vk(6>) = - f 7l ' 1 + 72 ' 2 ^ =U(6>) , (3.1) 

\72,1 -7i,2 / 

which was derived by Kaiser (1995) by combining appropriate combinations of third deriva- 
tives of the deflection potential ip and using (2.3). Its validity can also be checked by 
evaluating the gradient of (2.6b) and using that Aqi (— -^ln \0 — = 5(6' — 0). Using 
the definition (2.8) of g, Kaiser (1995) also obtained the relation 

VK{9) = - I J 1 + 91 / 2 )( Sl ' 1+ff2 ^u(0) , (3.2) 

l-<? 2 -2!V 92 1-01/ V 02,1 -01,2/ 

where 

K{6) := ln(l - «(0)) . (3.3) 

These relations show that if 7 is assumed to be known, n can be determined only up to an 
additive constant. However, only for weak lenses (k <^ 1) one can use the approximation 
7 g. From g, K can be determined up to an additive constant, or (1 — n) can be determined 
only up to a multiplicative constant. 

In Paper III, (3.1) was solved by integrating Vk over curves l\(t, 0) connecting a point 
b(A), < A < A, of the boundary of U with the point 0, and then averaging over all points 
of the boundary, 

K(0) = l^dA^^t^^.U(l A (t;0)) + i^dAK(b(A)) . (3.4) 

Since the second term is independent of 0, this equation yields k(0) up to an additive 
constant (i.e. , the mean of n on the boundary of U). Still there is much freedom left in the 
choice of the curves \\(t, 0). For a particular choice of these curves, results were presented 
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in Paper III. If the curves l\(t, 0) for fixed are chosen such that they do not intersect and 
that they cover the whole region U, then there is a one-to-one relation between (t, A) and a 
point 0' = l\(t, 0). In this case, we can rewrite (3.4) as 



k{0) 



dl A (*,0) dl A (i,0) 

X 



dA 



dt 



where we have defined for two 2-dimensional vectors a and b, 

a x b = ai&2 — G2&i 



dA «(b(A)) , (3.5) 



(3.6) 



In Paper III the curves 1 where chosen in such a way that the integrand in (3.5) agrees locally 
with that of (2.6b), which led to the fairly artifically looking integration contours. 

By generalizing the preceding method, KSFWB have noted that the starting points of 
the curves 1 need not be confined to the boundary of 14. One can proceed as follows: define 
curves l(t; 0, 0o) with 1(0; 0, 0o) = 0o, 1(1; 0, 0o) = 0, i- e -> which connect the points 0o and 
0. One can then integrate the differential equation (3.1) along each such curve. Averaging 
over all starting points, with an arbitrary weight function w(0q), one obtains 



k{0) 



J d 2 6 w(6 ) J 1 d* dl( ^ o) • U (l(f ; 0, O )) 



j u ^e oW (0 o ) 

+ / d 2 6 w(0o) k(0 o 
Ju 



(3.7) 



KSFWB have considered the special cases that w = 1 and that the curves l(i; 0, 0q) are 
straight lines (which implicitly assumes that the region U is convex). In this case, the 
constant, represented by the second term in (3.7), becomes the mean surface mass density 
R over the region U. Bartelmann (1995) has given an explicit equation for the resulting 
inversion equation, 

K (0)-R = [ d 2 6'U B (0',0)-V(0') , (3.8) 
Ju 

with 



1 



H B (0',0) = —- 



2A 



R 2 (0',0) 

\0-0'\ 2 



1 (0-0') 



(3.9) 



where R(0' , 0) is the length of a line segment from to the boundary of U passing through 

0'. 

All the preceding equations are also valid if we replace k(0) by K{0) - see Eq. (3.3) 
- and U(0) by u(0). In order to apply these inversion formulae, one needs to construct 
a continuous, different iable (complex) function g{0) from observed galaxy images, so that 
u(0) can be obtained from differentiation. In Sect. 5 we shall describe a simple method for 
obtaining such a smooth function. We also note that Eqs. (3.5) and (3.8) can be integrated 
by parts to remove the derivatives from 7; in the case of (3.8), this yields (2.11), with 
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V B (0',0) = n 



dH B dH B \ + . fdH B + dH B 



88[ 86' 2 



(3.10) 



due to the fact that H B (0',0) vanishes for 0' on the boundary of U. For a more general 
kernel vector H which does not vanish on the boundary of 14, this partial integration will 
yield boundary terms. Owing to the nonlinearity of u in g, this partial integration cannot 
be carried out in the equations for K. 

We can now easily check that V in (3.10) satisfies the constraint equation (2.12). In- 
serting (3.10) into (2.12) yields [with H B = H B (0',0)] 

Q(0,tf) = - / d 2 e'[(H B 1 -H B 2 )v 1 (0'-^) + (H B 2 + H B 1 )v 2 (0'-^)} 

71 Ju 

= — I d 2 8' [H B (V 1A + £> 2 , 2 ) + H B (2> 2il - V 1 , 2 )] , 
n Ju 

where we used again that H B vanishes on the boundary of U, and all indices separated 
by commas denote partial differentiation with respect to 6[. We then note that V can be 
written as second partial derivatives of the function f(0) := In |0|, as 



V 1 (0) = 

from which it follows that 
Z>i,i(0)+2?2, 2 (0) = 
Z>2,i(0)-2?i, 2 (0) = 



/,22(0)-/,ll(0) 



V 2 (0) = -f, 12 (0) , 



"2 (/.HI /.122W) 



- (Z\/) , = -7r^-<5(6>) , 



-\ (f, 112(0) +f, 222(0)) = ~\ W), 2 



-7T- 



5(0) , 



so that, after another partial integration, we obtain 

rB/„q a\ auB 



Q(0,0) 



dH*(ti,0) + dH»(0,0) 



01*2 



©(*) 



(3.11) 



Using the differentiation of H B as described in Bartelmann (1995), we finally find that 



Q(0,tf) = f or tf^e . 



(3.12) 



On the other hand, 



f d 2 d Q(0, •&) = -( d 2 d V • H B (tf, 0) = , 
Ju Ju 



(3.13) 



again because H vanishes on the boundary. Combining (3.12&3.13), we now see the validity 
of (2.12). 
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Using Eq. (3.7) to derive an inversion formula both has its strength and weakness, both 
essentially for the same reason: there is so much freedom left in the choice of the weight 
function w(0o) and the integration curves 0, 0o) that it will be difficult to single out 
'the best' choice for a given problem. However, this large freedom of choice may be used 
profitably, e.g., if the data set has 'holes' (which can be generated, for example, by some 
bright cluster galaxies which do not allow to measure g in a certain region of the cluster) to 
avoid such regions. We shall not discuss this issue further in this paper, but assume that 
the data for g are available on the entire region U. 

4 Finite field kernels: Noise filtering 

In this section we derive another inversion equation, again starting from (3.2) [or (3.1); we 
shall derive the inversion equation from (3.2), but can afterwards replace K by k and u 
by U to obtain an inversion equation of (3.1)]. The starting point of our considerations is 
the simple observation that the vector field u(0), obtained from observed distorted galaxy 
images, will in general not be a gradient field, but due to the noise, caused by the discreteness 
of galaxy images and their intrinsic ellipticity, it will also contain a rotational component. 
If u(0) is not a gradient field, equation (3.2) does not have a solution K(0) ! On the other 
hand, if u(0) is a gradient field, then all finite-field inversion equations are equivalent! [Of 
course, they may differ in their numerical applications; for example, the function (3.10) is 
not continuous if the region U has corners and must therefore be applied numerically with 
great care to avoid artificial features in the resulting ^-distribution]. If u is not a gradient 
field, the different finite-field inversion equations differ in their treatment of the rotational 
part of u. In other words, in this case the result of equations like (3.4) and (3.7) depends 
on the choice of the integration paths 1. 

Let us decompose the field u(0) into a gradient and a rotational part, 



where s(0) is a scalar field, and where we have defined the 'rotation' of a scalar in the 
second step, which is the gradient of the scalar field rotated by ir/2. Unfortunately, this 
decomposition is not unique, as shown by the simple example 



In order to determine VK from u, we thus need additional specifications for the decom- 
position. Since the rotational part of u is due to noise, it appears natural to perform the 
decomposition such that the rotational part is in some sense 'minimized'. In particular, we 
require that rots(0) should vanish if u is a gradient field. Also, we require that there is 
no systematic rotational component over 14, i.e., that the vector rots, averaged within 14, 



u(0) = VK{0) + rot s(0) = VK{0) + 



ds 

002 

-ds 



) 



(4.1) 
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should vanish. With these two conditions, the decomposition (4.1) is now uniquely deter- 
mined: from differentiation of (4.1), we obtain 

where A is the Laplace operator. If we require that 

s = on the boundary dU ofU, (4.26) 

then both of the above requirements are satisfied, and rot s is uniquely determined. We 
could of course choose any constant in (4.2b) without changing rot s. In fact, we shall not 
need to solve the boundary- value-problem (4.2), but only make use of the existence of its 
solution. 3 Then, we identify VK with X7K; although these two vector fields will be different 
in reality, since also the gradient part of u includes noise, we cannot separate this noise 
component from a true signal. The above identification is motivated by the expectation 
that the difference between VK and VK will be random, with zero mean over the field 14, 
and spatially uncorrelated on scales larger than the smoothing scale introduced later (see 
Sect. 5). Noting that K{6) can be determined only up to an additive constant, which we 
choose to be the average K of K over U, and that K and u are linearly related, we can 
make the ansatz 

K(0) -K = f d 2 6' H(0',0) -u(0') . (4.3) 
Ju 

We now replace u(0') by its decomposition (4.1) and integrate by parts; this yields 



K(0)-K=<f d\K(0')n(0') -H(0',0) + <f dA • H(0', 0) s(0') 

JdU JdU 

- f d 2 0' K(0') V • H(0', 0) - f d 2 6' s(0') V x H(0', 0) , 
Ju Ju 



(4.4) 



where dA is the length element on dU, and all differential operators operate on 0'. Here, 
n is the outward directed normal at the boundary of U over which the first two integrals 
are taken. We now consider the four terms in turn: the second term vanishes due to (4.2b). 
The last term can be made to vanish if we require H(0', 0) to be a gradient vector field with 
respect to 0', 

H(0',0) = V£(0',0) . (4.5) 
The first term can be made to vanish if we require 

H(0', 0) • n(0') = on the boundary of U. (4.6) 



It would be possible, of course, to solve (4.2) for s and to convolve the 'cleaned' data u — rot s with 
any finite-field kernel; this would yield the same resulting mass distribution as our new inversion 
formula derived below. However, one then has to repeat this 'cleaning procedure' for every new data 
set. Instead, it is our goal to develop a kernel which filters out automatically the uniquely specified 
rotational component and simultaneously yields a reconstructed density field. One then has to calculate 
the inversion kernel only once for a given geometry of the data field. 



11 



Finally, the third term can be made to agree with the left hand side of (4.4) if we require 




(4.7) 



where again A is the area of the region U. Eqs.(4.6) and (4.7) together define a Neumann 
boundary problem, which has a unique solution for C(0' , 0), up to an additive constant, i.e., 
it has a unique solution for H(0', 0). In the appendix, we present a closed- form expression for 
H for a circle, and give a fairly detailed description for the calculation of H for a rectangle. 

Finally, we want to point out that the vector field H needs to be calculated only once 
if similar field geometries are considered, due to the following scaling relation (which can 
be easily verified, using the explicit equations in the appendix) : if L denotes a scale of the 
data field (e.g., the side length of a rectangle or the radius of a circle), we introduce scaled 
coordinates x as = Lx, 9' = Lx.' . Then, 



5 Comparison of inversion methods 

In this section, we apply the various inversion techniques discussed in this paper and in Pa- 
per III to synthetic data. We choose a lensing mass distribution, distribute sources randomly 
in position, and assign an ellipticity to them according to an assumed intrinsic ellipticity 
distribution. For each such source, we then calculate the 'observed' ellipticity from the in- 
trinsic one and using the lens model. The data field to be analyzed thus consists of galaxy 
positions and ellipticities. Applying the four inversion techniques yields four reconstructed 
-^-distributions, which are compared to the true field, K true = In (1 — K true ). The difference 
K — K tvue is decomposed into its Fourier modes. We repeat this analysis for different source 
distributions many 50) times to obtain the power spectrum of the reconstruction-error 
field. Throughout this section, all angles are measured in arcminutes. The data field is 
assumed to be a square of side length L. 

5.1 Method 

5.1.1 Lens models 

The lensing mass distributions which we use here are constructed from superpositions of 
simple components of the form 




(4.8) 



where H is the kernel for L 



1. 




3/2 



(5.1) 
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where k,q is the central surface mass density, Q the center of the mass component, and C 
its core radius. Note that for \6\ ^> 9 C , the density behaves like an isothermal sphere. The 
Einstein radius of the corresponding singular isothermal sphere (i.e., that with the same 
behaviour for \9\ 3> 6 C ) is 6>e = Ko@ c - Four different mass distributions are investigated: 
Lens B is a strong, but sub-critical lens mass distribution consisting of two components of 
the form (5.1), e.g. a bimodal cluster or a cluster and a group of galaxies along the line of 
sight, and its center of mass is near the center of the data field; lens A is simular to B, but 
displaced to the boundary of the field; lens C is a weakly lensing mass distribution consisting 
of 3 components (e.g., groups of galaxies, or poor clusters) of the nearly isothermal mass 
distribution (5.1); two of the components are centered near the middle of the data field, one 
is at its left edge; lens D is the same as C but the third mass distribution at the left edge 
is removed. For each lens model the positions of the components, their core radii and their 
central surface mass densities are displayed in Table 1. To relate their lensing strength to 
that of a SIS-model we have added the corresponding velocity dipsersion of the SIS-lens in 
the table. For all models a surface- and contour plot of the field —K(x) = — In (1 — k(x)) is 
shown in Fig. 1. 

Table 1. Parameters of the four lenses A, B, C and D: Lens A, B and D are superpositions of two (a and 
(3) components of the form (5.1), lens C consists of three components (a, (3 and 7). The first three lines 
contain the peak positions Oq of the mass components, where the lower left corner of the field is the origin 
of reference, and all length scales are in units of arcminutes; the lines 4 to 6 (7 to 9) contain the core radii 
(central mass densities) of the superposed mass distributions; in the last three lines the velocity dispersions 
of the corresponding isothermal sphere models are indicated. 
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Fig. 1. Contour- and surface plot for the function —K(x) = — In [1 — t(x)], where k is the two-dimensional 
surface mass density of the lens A,B,C and D respectively. In each panel, the size of the field is 7. '5 X 7. '5, 
K is calculated on a 50 X 50 grid; the spacing in the contour lines is 0.1 for lens A and lens B and 
0.05 for lens C and D 



5.1.2 Generation of the source data and calculation of the 'observed' data field 

We distribute sources randomly in 6 (i.e. on the observer and not on the source sky), which 
is a valid procedure due to the fact that the magnification bias basically vanishes for the 
faint blue galaxy population, since their source counts behave nearly as N(> S) oc 
Their density is no/(arcmin) 2 , and typically no = 50. The intrinsic ellipticity e' s ' of the 
sources is defined as in Paper III. [Briefly, e' s ' is a complex number whose modulus is 
(1 — r)/(l + r), for an elliptical source with axial ratio r < 1; for general brightness profiles 
of the sources, e' s ' is defined in terms of the tensor of second brightness moments. The 
phase of e' s ' describes the orientation of the source.] In our simulations, the intrinsic 
ellipticity e' s ' was drawn from a Gaussian distribution of the form 

Ps(€ (s) )= m 1 - 1/p2 , e-' e(S) l V > (5-2a) 
7T/CT (1 — e l /p ) 

where p s (e^) d 2 e' s ' is the probability that the source ellipticity lies within d 2 e' s ' of e' s '. 
Hence, the quantity p controls the width of the intrinsic ellipticity distribution, and we 
expect that with increasing p, the reconstructed mass density will become noisier. 

Other authors have used ellipticity distributions different from (5.2a); in Papers IMI, 
and in KS, Gaussian distributions in the absolute value of 

1 + |e| 

were used (see Paper II, eq. (2.22) and KS, page 445, second column, and note that the 
quantity \e\ used by KS is equal to |%| defined in Paper II). Hence, these authors used an 
ellipticity distribution of the following form: 

»fr' ,) >= ,i P (l-e-V») «- W ' ^ 

Using that for sources with elliptical isophotes |e' s '| = and |x' s '| = ^~^ 2 , we 

have calculated the probabilities p p (r) and pn(r) that a source drawn from an elliptic- 
ity distribution (5.2a) or (5.2c) has an axis ratio r; the widths R and p were chosen as 
R, p e [0.1, 0.2, 0.3, 0.4, 0.5]. The results can be seen in Fig. 2: the solid curves, p p (r), were 
obtained from (5.2a), where with increasing p the maximum is shifted to the left, i.e., 
elliptical sources become more likely; the analogous curves, pn(r) for (5.2c) are dashed. 
One can see from Fig. 2 that for p <; 0.2, p p {r) ~ PR=2 P {r); this is due to the fact that 
for small e, (5.2b) becomes % Rj 2e. In particular, this implies that the simulations in KS 
and in Papers I&II, where an ellipticity distribution of the form (5.2c) with R = 0.3 was 
used, consider only rather circular sources; this would correspond to a width of p = 0.15 
in our simulations. 

Finally, the 'observed' ellipticity for each galaxy is calculated from the lens model, 
using 

e = fr- . 5.3 

We want to stress that this transformation law is valid only for non-critical clusters, where 
\g\ < 1 everywhere. In the more general case of critical clusters, one should use the 
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0.2 0.4 0.6 0.8 1 

r 



Fig. 2. Solid lines: probabilities p p (r) that a galaxy drawn from the ellipticity distribution (5.2a) has an 
axis ratio r, where p £ {0.1,0.2,0.3,0.4,0.5}. Dashed lines: probabilities prt(r) that a galaxy drawn from 
the ellipticity distribution (5.2c) has an axis ratio r, where R £ {0.1,0.2,0.3,0.4,0.5}. For p < 0.2, the 
probabilities p p (r) and PR=2p( r ) roughly coincide. 



ellipticity parameter %, as defined in Paper I. Thus, the synthetic data set consists of a set 
of values for the galaxy positions g i and their corresponding ellipticity ej. 

5.1.3 Determination of g, 'inner-smoothing' 

If we average the observed ellipticities of the images over a region of the lens within which 
g is approximately constant, their expectation value is 

(e) = -9 , (5.4) 

due to the fact that the sources are intrinsically randomly oriented. This fact has been 
shown by Schramm & Kayser (1994) and can be obtained simply from an angular integra- 
tion of (5.3). On a regular grid of points 

^• = (m,ja) , 0<i,j<N , (5.5) 

with a = L/N being the size of a gridcell, and N the number of gridpoints per dimension 
on the square, we have calculated an estimate of 9(0'-) = gij by a weighted average, 
according to (5.4), 

g l3 = - ^ Wmem , (5.6) 

with Gaussian weights 
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Ad? 



w m oc exp ( 13 A gm ] , (5.7) 



13 

and the smoothing scale A9 can be chosen appropriately. As described in detail in Paper II, 
it is useful to adapt the smoothing scale to the strength of the signal; i.e., in a region where 
\g\ is large, the ellipticity of the observed galaxy images will be dominated by the shear 
effect of the lens, and an accurate value of g can be obtained by averaging over a few 
images only, whereas in regions of weaker shear, the average should extend over a larger 
number of images. We take the same smoothing procedure as in Paper II, i.e., we choose 



AOij = A0 



o^-l^.il 2 )^ , (5.8) 



where the prefactor AQq should be of the order of several times the mean separation of 
galaxy images (e.g., A6 ?a -J=p), and the exponent (3 determines how the smoothing 

length decreases with increasing \g\, i.e., increasing \S\. For the inversions shown here, 
we chose (3 = 2. Since we later will also smooth the reconstructed density field, we call 
the smoothing of the ellipticities to get an estimate of g the 'inner smoothing' and that 
of the density field 'outer smoothing'. We point out that we applied the same smoothing 
procedure (with the same smoothing lengths) to all four inversion methods. Changing the 
smoothing lengths affects the quality of the reconstruction differently for each method, 
and each method will have a different optimal smoothing length with respect to a given 
quality criterium for the reconstruction. For instance, if the quality criterium is to resolve 
the heights of density peaks, then the maximally acceptable smoothing length of the KS 
inversion will be smaller than that of the other three inversion techniques. However, 
we want to emphasize that we do not develop an optimized inversion with respect to 
smoothing for each inversion kernel in this Paper, but we compare the noise sensitivities 
and systematics for the four kernels; in particular we show that the kernel developed in 
Sect. 4 is the least noise-sensitive finite-field kernel, and the increase in noise compared to 
the KS reconstruction is compensated by far by the loss of the boundary artefacts. 



5.2 Discretized form of the inversion formulae 



5.2.1 KS-inversion 

With the values of gij obtained according to (5.6), we calculate the KS estimate [see 
Eq. (2.10a)] of n iteratively, as suggested in Paper II 

2 N 

K n+1 (0'ki) = W * W J i 1 - M V * K - O'ij) 9ij] , n > . (5.9) 

i,j=0 

The starting value of this iteration is k° = 0, and convergence is obtained after only a few 
iterations. Here, Wi = 1/2 for i = and for i = N, and Wi = 1 otherwise. The resulting 
solution is then interpolated on the grid 

Bu = ( [A; - 1/2] a , [I - 1/2] a ) (5.10) 

by 
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K kl = K(0 kl ) = \ [K(0' kl ) + K(0' fc _l,l) + «(0'fc,I-l) + • (5-11) 

Since we use the nonlinear form of the reconstruction, i.e., we take g as the observable, not 
the shear 7, the result of the iteration determines K, as defined in (3.3), up to an additive 
constant. Hence, our estimate of K k i = K(0 k {) from the KS method, uncorrected for the 
finite field, is 

-K = \n(l-n kl ) , (5.12) 
where the constant K cannot be determined. 

5.2.2 Finite field kernels 

The other three reconstruction methods we want to investigate here are all based on the 
vector field u(0), defined in (3.2). We have calculated this vector field on the grid 6'- 
from (5.6), using finite differencing, thus obtaining the values u^- = u(0^-). Using (3.8), 
with k replaced by K and U replaced by u, we obtain a second estimate for K from the 
method developed by KSFWB and Bartelmann (1995), which we shall call for notational 
simplicity 

N 

Kf d -K = a 2 Y J WiW^iO'^Ou) ■ u {j , (5.13) 

i,j=0 

with H B given in (3.9). Similarly, the estimate from our kernel developed in this paper 
becomes, according to (4.3), 

N 

Kff -K = a 2 Y^ WiWjH^eu) ■ u^- , (5.14) 

i,j=0 

with the kernel H determined numerically, as described in the appendix. Finally, we obtain 
an estimate of K k i from the inversion method described in Paper III, which we shall term 

^ki- 

5.2.3 'Outer smoothing' 

Obviously, the smoothing introduced by AO is fairly moderate, if A6q is of the order of the 
mean galaxy separation, and the resulting estimates for K will be very noisy for reasonable 
values of the parameter p in the ellipticity distribution (5.2). We remind the reader that 
we needed the first smoothing to obtain a smooth function g(0) which can be differentiated 
(of course, g will be fairly noisy too). Hence, we apply a second smoothing on the resulting 
distribution of K, i.e., we define for all estimates the smoothed fields 

E^exp(H^-0 fc dV€n) 
K k i = t± ^ , (5.15) 



E^exp^i-^lVen) 



with the same smoothing length for the KS, SS, S and B reconstructions. In principle, the 
smoothing length sm can depend on the position k i, and an 'optimized' reconstruction 
could adjust this local smoothing scale by some appropriate measure, such as a local 
measure of a x 2 -statistics [see Eq. (4.8) of Paper II]. We shall not develop such a scheme 
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here. Instead, we will use 6 sm to be a constant, and to be typically of the order of the 
gridsize or larger (this removes grid effects). 

5.2.4 Power spectrum of the error-field 

We then compare the resulting smoothed field K from the various reconstruction methods 
with the true field K ty:ue = l n (l — ^true)- For this purpose, we define the Fourier transform 
of the difference between the reconstructed field K and the true field K tvue , 

£>(k) = ± f d 2 6 e ik [K{0) - K true (0)] , (5.16) 

where K stands for K KS , K ss , K s and K B , respectively. We then define the power 
spectrum of this difference to be 

P(k):=(D(k)D*(k)) , (5.17) 

where the average extends over the directions of the (appropriately binned) k- vectors and 
various realizations of the source distribution. 

5.3 Results 

5.3.1 Results for the power spectra 

For a given lens model, the power spectrum of the error field depends on several parameters: 
the number N of gridpoints, the size L of the data field, the density no of galaxies, the width 
p of the ellipticity distribution, the inner smoothing length 9q and the outer smoothing 
length 9 sm . For all power spectra shown here the number Z max of simulations which have 
been made to calculate the power spectrum is equal to 50 and N = 50. If not stated 
otherwise, L = 7.5, which is the size for a currently available CCD at the CFHT (Fahlman 
et al. 1994, Tyson 1994). The density of galaxies varies between 40 and 80 (per square arc 
minute), and 0.2 < p < 0.4 4 ; the inner smoothing length varies from 0.2 to 0.5 (arcminutes) 
and the 9 sm is typically of the order a/y/2. In order to emphasize the differences of the 
reconstructions at small length scales, we have plotted P(k)k 2 instead of P(k). 

We first consider reconstructions of the lens B, which is an example for a fairly strong 
but non-critical cluster. Since the ellipticity distribution of the faint background popu- 
lation is currently not well determined, we consider ellipticity distributions with p in the 
range of 0.2 to 0.4. For no = 50 and p = 0.2 we reconstructed the lens with 0$ = 0.25, 
#sm = (reconstruction Bl) and 6o = 0.34, 6 sm = a/y/2 = 0.11 (reconstruction B2); the 
corresponding power spectra are shown in Fig. 3a and 3b. 

One sees that the KS power spectrum has a large spike at k = 2n/L, which is due 
to the systematic artefacts near the boundary of the CCD field. Despite this, at larger 
wave numbers the power spectra of the KS, SS and S error look quite simular, where 

4 Recent investigations by Brainerd et al. (1995; Brainerd 1995, private communication) indicate that 
the observed (seeing convolved) galaxy images down to r ^ 26 are fairly circular: the probability 
distribution peaks at r = 0.88, whereas the behaviour for small r seems to be similar to the case of 
a probability distribution according to (5.2a) with a width of p = 0.2. Therefore, simulations with 
p = 0.2 seem to be quite realistic. However, since the ellipticity distribution for galaxies in the blue 
light is probably broader, we also present test simulations with larger values of p below. 
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Fig. 3. Power spectra of the error field of the reconstructed lens B, for parameters TV = 50, L = 7.5, 
no = 50, p = 0.2; the inner and outer smoothing lengths were varied: (a) 6q = 0.25, # S m = 0, (b) 6q = 0.34 
and 6 sm = 0.11. The solid (long-dashed, short-dashed and dotted) lines correspond to the KS- (SS-, S- 
and B-inversion) respectively. The spike in the power spectrum of the KS-inversion at k = ^ comes from 
the boundary artefacts of this reconstruction method; for larger values of k the power spectra of the KS-, 
SS- and S-reconstruction roughly coincide, the B-reconstruction is fairly noisy. For wave vectors k > 3, 
and for all inversion techniques, the power spectrum mainly reflects 'real noise' in the inverted i^-field; in 
other words, for a small smoothing length only a minor part of the power spectrum at large k is due to 
systematic effects, such as not resolving small scale maxima of the density field due to smoothing. The 
exponential cutoff at k ^ 8 in Fig. 3b is caused by the 'outer smoothing', sm = 0.11 



the quality of reconstruction is slightly decreasing from the KS- to SS- and S-method 
at those wave vectors; the power spectrum of the B-inversion is considerably worse, i.e., 
the reconstruction is rather noisy on small scales. The reconstruction Bl was obtained 
without any 'outer smoothing', 9 sm = 0; the influence of an outer smoothing can be seen 
by comparison with Fig. 3b. In this model, the inner smoothing length was also larger, 
which yields a general decrease of all power spectra. The outer smoothing becomes visible 
for k J> 8 by an exponential suppression of the power spectra. For the reconstruction B2 
we also show a single realization of the KS-, SS-, S- and B-inversion in Fig. 4: we have 
used the same source distribution for all inversions and all K-estimates are shifted such 
that the average over the field equals the average of K true . 

Inspection by eye indicates that the KS-reconstruction is the smoothest one, and 
thus shows the least amount of noise; one can see as well that the SS- and S-reconstruction 
have roughly the same quality, and that the B-reconstruction is by far the noisiest. For the 
smoothing length used, the KS-inversion does not resolve the maxima as well as the three 
finite-field kernel inversions do. The KS-reconstruction shows the well known artefacts 
at the boundary, i.e., local maxima at the corners of the field as well as local minima at 
the middle of the edges. The other reconstruction methods also show enhanced errors at 
the boundary; however, they are not systematic, but caused by increasing noise, since for 
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Fig. 4. Surface and contour plots of — J4T KS , — K ss , —K s and — obtained by the inversion of the lens 
B using the parameters of reconstruction B2, i.e. L = 7.5, N = 50, no = 50, p = 0.2, 6q = 0.34 and 
0sm = 0.11; the contours differ by 0.1 21 



positions near the CCD edge, the quantity g is determined by averaging over fewer galaxies 
(in the extreme case, at a corner, one averages only over one fourth of the galaxy number 
that one averages for a point away from the boundary). In the reconstructions shown in 
Fig. 4, one can also see the impact of intrinsically very elliptical galaxies, see, e.g., the 
right front corner of the 'CCD'. Near that corner there is an intrinsically very elliptical 
source yielding a large estimate for the distortion; therefore, all reconstructions show a 
local maximum at this corner. For the KS-reconstruction, this local maximum is not so 
pronounced, because the reconstructed field is already strongly biased at the corners, and 
it appears that the KS-kernel is less sensitive to individual intrinsically elliptic sources 
than the three finite-field kernels. 

On second sight, however, the impression that the KS-reconstruction is less noisy 
than the SS- and S-reconstruction must be modified. Namely, this impression is mainly 
due to the reconstruction near the boundary. As we stated before, due to the fact that 
we have less information about the distortion near the boundary, compared to points 
closer to the center of the field, we expect that all finite field inversions will be noisier 
close to the boundary. For the KS-reconstruction, the reduced information about the 
distortion near the boundary is substituted by artificially setting the shear to zero outside 
U. Hence, one uses more (but wrong) information in the KS-inversion, which yields a 
smoother appearance of the reconstructed density field. 

We next study the reconstruction with the four techniques for an ellipticity distri- 
bution with a width of p = 0.3, where more than 50 percent of the sources have an axis 
ratio smaller than 0.6. We have reconstructed the lens B for a galaxy density of no = 50, 
ellipticity width p = 0.3 and outer smoothing length of sm = 0.11, and we varied the inner 
smoothing length: # = 0.25 (reconstruction B3), 9o = 0.3 (reconstruction B4), # = 0.33 
(reconstruction B5), 0$ = 0.51 (reconstruction B6). The resulting power spectra can be 
seen in Fig. 5. 

As expected, the power spectrum becomes larger if one compares a reconstruction 
for an ellipticity distribution with width p = 0.2 with that of a width of p = 0.3 using 
a fixed smoothing length, see, e.g., Fig. 3b and Fig. 5c: To first order, the power spec- 
trum is roughly doubled for all inversion kernels; however, comparing Fig. 5c and Fig. 3b 
in more detail shows that the sensitivity to noise, introduced by the more elliptical source 
distribution, increases from the KS- to the SS- and S-reconstruction, and it is strongest 
in the B-reconstruction. Fig. 5 also shows how the power spectrum changes with increas- 
ing smoothing length, with all other parameters kept fixed. Generally, one needs a larger 
smoothing length for the SS-, S- and B-inversion than for the KS-inversion to lower the 
noise, where always the SS-kernel yields the best reconstruction of the three finite-field 
kernels. The short wave vector power spectrum for KS is increasing with larger smoothing 
length, since for KS this smoothing length already lowers the two maxima of the recon- 
structed if-field considerably. For the finite-field kernels this suppression of density peaks 
becomes important at larger smoothing length than for KS, presumably because for the 
KS-inversion one has to use the quantity g itself, whereas for the finite-field kernels one 
needs the derivatives of g. Finally, using a very large smoothing length like 9q = 0.5 or even 
larger yields a convergence of the power spectra of the 4 inversion types for large wave 
vectors: then, the reconstructions become all rather smooth, but small scale structures 
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Fig. 5. Power spectra of the error field for reconstructing the lens B with no = 50, p = 0.3, # S m = 0.11; 
the 'inner smoothing' length was: (a) 6 = 0.25, (b) 6 = 0.3, (c) 6 = 0.33 and (d) 6 = 0.51. In all 4 
panels, the solid (long-dashed, short-dashed and dotted) line corresponds to the KS-power spectrum (SS-, 
S- and B-power spectrum). As expected, the smoothing length has to be increased compared to the case 
of p = 0.2 to get a noise-poor reconstruction, i.e. a small power spectrum. 



can no longer be reconstructed successfully. In Fig. 6 we demonstrate how the quality of 
reconstruction changes with decreasing and increasing galaxy density. 

Reducing the galaxy density from no = 50 to no = 40 means that the mean separa- 
tion increases by about ten percent; hence, keeping the smoothing length constant (i.e., 
averaging over fewer galaxies to estimate g) yields an increase in the power spectrum for 
all inversion methods, as can be verified by comparison of Fig. 5b with Fig. 6a. However, 
the loss of reconstruction quality increases from the KS-reconstruction to the SS-, S- and 
B-reconstruction. An increase of the smoothing length to 0.4 and 0.5 can reduce the power 
spectrum again. Clearly, increasing the galaxy density reduces the power spectra, and the 
reconstruction with the KS-, SS- and S- technique become very simular, where, of course, 
the latter two do not show the spike of the power spectrum at small k; thus, at galaxy 
densities of 80 the SS- and S-reconstruction techniques are as good or superior to the 
KS-reconstruction on all scales. 

If the ellipticity distribution of the sources is more extreme, p = 0.4, then for a galaxy 
density of 50, the three finite-field kernels yield only very noisy reconstructions, compare 
Fig. 5a and 7a. Using a large smoothing length of 6 sm = 0.5, see Fig. 7b, improves the 
reconstruction considerably and yields similar results as for the KS-reconstruction, but at 
the cost of small-scale resolution. The first maximum in the power spectra of the three 
finite-field kernels in Fig. 7b comes from the fact that with very elliptical sources, the 
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Fig. 6. Power spectra for the error field where the lens B was reconstructed, using a width of p = 0.3 in 
ellipticity distribution and # sm = 0.11; galaxy density and inner smoothing lenght were chosen as follows: 
(a) n = 40, O = 0.3, (b) n = 40, O = 0.4 (c) n = 40, 6 = 0.5 and (d) n = 80, 6 = 0.3. Reducing the 
galaxy density and not increasing the inner smoothing length yields a increase of the power spectrum for 
all 4 inversion techniques, compare Fig. 5b and Fig. 6a; the loss in quality of the reconstruction increases 
from the KS-, to the SS-, S- and B-reconstruction 



reconstruction errors near the edges and especially at the corners become very large. In 
contrast to the KS-inversion, where the reconstruction at the outer parts of the CCD is 
systematically affected, one obtains for the finite-field kernels a 'wavy density field' at the 
outer parts. Fig. 7c demonstrates how much an increase in the galaxy density helps in 
the improvement of reconstruction quality if the galaxies are intrinsically very elliptical: 
compare 7a and 7c, where the parameters for which the power spectra are calculated differ 
only in the galaxy density. Also, as expected, if the galaxy density is large, e.g., no = 80, 
a change in the ellipticity distribution towards more elliptical sources does not affect the 
reconstruction as much as for the case of a lower galaxy density (see the change in the 
power spectra from Fig. 6d to 7c, where for no = 80, both p and 6>o are increased from 
0.3 to 0.4, and compare this to the change of the power spectra from Fig. 5b to 7a, where 
the galaxy density is 50, and also both p and 6q are increased from 0.3 to 0.4). We want 
to point out that the loss in reconstruction quality for the finite-field kernels for p = 0.4 
compared to p = 0.3 is mainly caused by the increased fraction of galaxies which have 
an intrinsic axis ratio of 0.2 or smaller (see Fig. 2). However, since spiral galaxies have a 
finite disc, even an edge-on galaxy will have an axis ratio of 0.2 or larger; therefore, an 
ellipticity distribution with a width of p = 0.4 appears quite unrealistic. Despite of this, 
Fig. 7 shows that cluster inversion can be done successfully even in the presence of very 
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Fig. 7. Power spectra for the extreme case of an ellipticity distribution with a width of p = 0.4; this 
ellipticitiy distribution is not realistic since quite many galaxies with an axis ratio of 0.2 or smaller are 
generated; however, the power spectra demonstrate that the KS-inversion method can extract the distortion 
quite well even in the presence of extremely elliptical galaxies. The galaxy density and inner smoothing 
length have been varied as follows: (a) no = 50, 6q = 0.4, (b) no = 50, 6q = 0.5, (c) no = 80, 6q = 0.4 and 
(d) no = 80, 8q = 0.5; in all four panels, the outer smoothing length was fixed as # sm = 0.11. As in the 
preceding figures, solid (long-dashed, short-dashed and dotted) lines denote KS- (SS-, S- and B-) power 
spectra 

elliptical sources, if only the observed galaxy density is sufficiently high. 

To see what happens if the inverted mass distribution is not centered on the center of 
the CCD, we reconstruct the mass distribution of lens A. This does not mean that we expect 
an observer to position the telescope on purpose in such a way; however, it demonstrates 
how large the systematic error of the KS-reconstruction can become in unfavourable cases; 
in addition, we note that in some cases (see Bonnet, Mellier & Fort 1994) where bright stars 
are in the foreground, one sometimes must position the CCD in a non-optimal position for 
cluster inversion. 

Fig. 8 demonstrates that the KS-reconstruction quality decreases very strongly, and 
that for the finite-field kernels this reduces the reconstruction quality only slightly (due 
to our smoothing procedure, the noise of the reconstruction at the boundary is larger). 
A more realistic observational situation is simulated in the reconstruction of lens C and 
D. It may happen that someone tries to invert a mass profile of two groups of galaxies 
(lens D); then, a nearby additional dark mass distribution (lens C) will lead to large errors 
in the KS-reconstruction. (Such dark mass distributions do exist as has been proven 
by the detection of shear fields around bright 1-Jansky quasars, caused by 'dark' mass 
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Fig. 8. Power spectra for the error field for the reconstruction of the lens A: in all panels, no = 50 
and # sm = 0.11; the ellipticity width and inner smoothing length are varied. Since the mass profiles of 
lens models B and A are quite similar, these power spectra are shown to demonstrate the increase of 
reconstruction error for the case where the center of the mass distribution is shifted to an edge of the CCD 
field. Hence, compare Fig. 8a with 3a, 8b with 5b, 8c with 7a and 8d with 7b 



distributions; B. Fort, private communication.) 

For lens D, the KS-reconstruction yields a good estimate, since the shear outside the 
data field, which is set to zero in the KS-inversion, then is indeed very small. Hence, the 
spike at small wave vectors in Fig. 9a is very small, and the KS-, SS- and S-estimates of the 
mass distribution are then almost equally good (Fig. 9a) . The additional mass distribution 
at the boundary of the CCD in lens C, however, yields again a prominent rise of the spike 
(Fig. 9b). Hence, also for the reconstruction of a very weak cluster or groups of galaxies, 
it is profitable to use a finite-field kernel inversion. 

In order to demonstrate at which positions the errors in the 4 different inversion 
methods are particulary large, we have plotted the root mean square of the reconstruction 
error as a function of position, see Fig. 10. 

The rms error was obtained by 500 simulations of reconstruction of lens B, using a 
galaxy density of 50, a width in the ellipticity distribution of p = 0.3, an inner and outer 
smoothing length of 0.3 and 0.11. The corresponding power spectra of the error fields are 
those in Fig. 5b. The rms error at a position is defined as follows 

AK(t) :=^< (K shiit (0) - K tIue (0)) 2 > , 
where -PC s hift is the reconstructed K-field obtained for the estimator considered, shifted 
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Fig. 9. Power spectra of the error field for reconstructedlens D (a) and lens C (b); as in previous plots, the 
solid line (long-dashed, short-dashed and dotted lines) belong to the KS- (SS-, S- and B-) power spectrum, 
(a) For a weak mass distribution at the center of the CCD field, the KS-reconstruction yields a good 
estimate, where the finite-field kernel reconstructions are a little bit noisier on small scales; (b) if there is 
an additional mass at the boundary of the CCD, the quality of reconstruction decrases considerably for 
the KS-reconstruction, whereas it stays approximately the same for the finite-field kernel reconstructions 

such that the sum of -ftT s hift over the data-field equals that of K trU e- The rms error field 
for the KS-reconstruction is shown in the upper left panel of Fig. 10: i) the reconstruction 
error is largest at the corners where it exeeds 0.2; ii) the rms-error field shows a maximum 
at the center of the data field; iii) 'behind' the central maximum, there is a strong increase 
of the error towards the edge of the data-field. The reconstruction error has to be large at 
the corners and it is large along those edges where the shear field is still high (this explains 
why the rms error is not as large on the remaining three edges). However, we want to 
point out that the comparison done in Fig. 10 is not 'fair', due to the normalization of 
the K-field. Since the normalization is affected by the boundary artefacts, it may increase 
the error in the middle of the data field. However, we have not come up with a 'fairer' 
normalization. For this reason, the comparison of Fig. 10a with the other panels should be 
considered with care. 

The finite-field kernels rms error field are qualitatively similar to each other: at the 
edges and corners the rms error is also increased, but contrary to the KS-reconstruction, 
the rms error arises only from the increased noise (smoothing over fewer galaxies) in the 
mass distribution at these positions. This can partly be avoided by a more elaborated 
smoothing procedure, which, however, will not be discussed in this paper. 

A comparison of the rms errors of the finite-field kernels with each other yields the 
same ranking for these inversion techniques as from considering the power spectra: the 
size of the rms error increases from SS-, to the S- and B-reconstruction. In particular, the 
error of the SS-inversion seems to be spatially flatter than for the other two finite-field 
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Fig. 10. Root mean square error of the reconstruction as a function of position for the KS-, SS-, S- and 
B-reconstruction. The number of simulations is 500, no = 50, p = 0.3, 6q = 0.3 etna C7sni 

= 0.11; the 

g^responding power spectra for the error field are displayed in Fig. 5b. The vertical axis is drawn from 
to .2; the contour lines are in the same interval, with a step size of 0.0125; the size of the CCD field is 7.5 
arcminutes, as usual 



inversion techniques. 
6 Discussion 

In this paper we have derived a cluster inversion formula which yields the surface mass 
density distribution from observations of the image distortion on a finite data (CCD) 
field. This new inversion equation is obtained by convolving the vector field u(0), which 
is obtained from the observable g(0) and its derivatives - see Eq. (3.2) - with a kernel 
H (#',#), which in turn is obtained by solving a Laplace-like equation with Neumann 
boundary conditions. We have explicitly constructed the kernel H for a rectangular field, 
as described in the appendix (for a circular field, the kernel H is known in closed form; 
for more general field geometries, H has to be calculated numerically). We then have 
compared this new formula with the nonlinear version of the KS inversion equation and 
two other finite-field inversion formulae derived in Paper III, and KSFWB and Bartelmann 
(1995), respectively, by applying these various methods to synthetic data and performing 
a spectral analysis of the deviations of the resulting mass profile from the original mass 
distribution. The main finding of this comparison is that our new formula in all cases 
considered is better than the two other finite-field kernels, and at worst only slightly 
more noisy than the KS reconstruction; however, this slight increase in noise is more than 
compensated by the removal of boundary artefacts with which the KS method is burdened. 
If the mass distribution extends near to the field boundaries, then our new method yields 
reconstructions which deviate less, even on small scales, from the input mass distribution 
than that obtained from the KS formula, since then the latter can be seriously hampered 
by boundary effects. This effect becomes even more pronounced if rectangular data fields 
with side ratios not close to one are considered (which we have not done in this paper, 
but which was demonstrated in Paper III). In most cases, our new reconstruction kernel 
yields only slightly better results than the method of Paper III, whereas in nearly all cases 
it is considerably better than that of KSFWB and Bartelmann (1995). We also want to 
point out that the noise of our new reconstruction method seems to be uniform over the 
data field, except that it is slightly larger at the boundary (which must be the case since 
there is less information about the distortion near the boundary), whereas the other three 
inversions have their noise more concentrated towards the boundary. 

The finite-field inversions S and B were derived from line integrations of the gradient 
of the surface mass density (3.1) or (3.2) and averaging over many different curves 1; hence, 
they are special cases of the more general inversion formula (3.7). This equation leaves a 
lot of freedom (or rather: arbitrariness), i.e., the choice of the weight function w(0o) and 
the choice of the curves 1 connecting 0q an d 0. In Paper III, these integration contours 
were chosen such that the resulting inversion equations agrees locally with that of the KS 
formula; in KSFWB and Bartelmann (1995), this is not the case. As argued in Paper III, 
this requirement appears reaonable due to the singularity of the kernel for 0' close to 0, 
and we attribute the considerably worse noise properties of the B inversion relative to the 
S inversion to this fact. We also note that we have used the B inversion in the form (3.8), 
not in its integrated form; some test calculations have shown that the latter one yields 
results which are noisier than those obtained from (3.8), due to the unsmooth behaviour of 
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R(0' , 0) - see Eq. (3.9) - in the direction towards the corners. It may well be that for field 
geometries with smooth boundary curves this difference will vanish; however, we think 
that the local anisotropy of the kernel H B for 0' —> is the main cause for its increased 
noise relative to the S inversion. The choice of the integration curves for the S inversion 
in Paper III appears artifical at first sight, but has proven to be reasonable. However, for 
both the S and the B inversion (in their present form) it is essential that the data field is 
convex. 

The rationall behind our new inversion method was the consideration that the 'ob- 
served' vector field u(0) will in general not be rotation-free. Our new kernel was derived by 
requiring that any rotational component of u is filtered out, and that the decomposition of 
the field u into a gradient and a rotational part is such that the latter is minimized in a par- 
ticular sense (see Sect. 4). These two requirements then specify the integral kernel H(0' , 0) 
uniquely. In particular, the generalization to a different field geometry is reduced to the 
solution of a Laplace equation with Neumann boundary conditions, and non-convex data 
fields can be treated with the same method. We shall treat other geometries elsewhere, in 
particular that of the (non-convex) WFPC-2 field on board HST. 

The fact that our new inversion method has not significantly worse noise properties 
than KS is surprising at first sight, given the fact that it needs the data in a differentiated 
form. Of course, the field u is more noisy that the input field g; however, the integration 
(4.3) appears to cure this additional noise caused by finite differencing. 

We have not tried to find optimal smoothing procedures in this paper, and in fact 
used the same smoothing procedures for all four inversion techniques. It is likely that each 
of these inversions has its own optimal smoothing scale, but that will depend strongly on 
the quality of the data available. Hence, the choice of the smoothing procedure has to be 
made after the actual data set is obtained. 

Looking at the power spectra in Sect. 5 one might suggest alternative inversion meth- 
ods by combining our new formula with that of KS. This could be done either by a Fourier 
decomposition of the KS reconstruction and the SS reconstruction, and then using the SS 
components for the large scales and the KS components for small scales. Alternatively, 
one could use the SS reconstructed density field to calculate its shear outside the data field 
and then apply the KS inversion, using the observed field g inside the data field and the 
calculated field g outside the data field. Several more variants of this kind can be thought 
of. On the other hand, it remains to be seen how useful the method of Broadhurst, Tay- 
lor & Peacock (1995) and, in particular, that of Bartelmann & Narayan (1995) will be. 
Certainly, both of these new methods (which make use of the magnification effect which is 
completely neglected in the inversion techniques discussed here) will yield additional infor- 
mation about the surface mass density of the cluster; in particular, the additive constant in 
the /^-distribution can probably be determined from these methods. It will be difficult to 
incorporate this additional information into inversion formulae of the type discussed here, 
and one will have to think about more general inversion techniques which can make use of 
all available information. Most likely, this will lead to some kind of maximum-likelihood 
approach for a parametrized mass distribution. Such a method should also take into ac- 
count the redshift distribution of the sources, and in fact should be able to determine that 
distribution if several clusters at different redshifts are simultaneously fitted. A very first 
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step in this direction was taken in Smail, Ellis & Fitchett (1995). 
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solve our Neumann problem, as done in the appendix, and M. Bartelmann and J. Ehlers 
for their most detailed comments on the manuscript. This work was supported in part by 
the Sonderforschungsbereich SFB 375-95 of the Deutsche Forschungsgemeinschaft. 



Appendix 

Here we want to derive the solution for the 'Poisson' equation 

AC(0',0) = -5(0-0') + j , (Al) 

for inside the region U with area A, together with the boundary condition 

H(0',0).n(0') = -± 7 £(0>,0) = O (A2) 
dn' 

for all points 0' on the boundary dU of U. Here, 

H(0',0) = V£(0',0) , (-43) 

all differential operators in this appendix are to be taken with respect to 0', and n is the 
outward-directed normal vector on the boundary of 14. The existence of a solution of this 
Neumann problem is guaranteed by the compatibility condition 

J d 2 0' AC{0',0) = , (A4) 

(see, e.g., Garabedian 1986, p. 230), and it is unique up to an additive constant, which 
implies that the solution for H is unique. Here we shall derive the solution for two special 
geometries, a circle of radius R and a rectangle (which is the most relevant case due to the 
shape of CCDs and most, though not all, CCD arrays). 

For a circle of radius R, the solution of the Neumann problem can be found in text- 
books (e.g., Garabedian 1986, p. 248); we only quote the result for H: 

„, 1/0' 0-0' R 2 0/\0\ 2 -0' . 

I |0-0 | R20/\0\ 2 -0' 

The validity of this solution can be easily checked. 

Consider next a rectangle with sidelengths L\ and L2. In this case, the solution of 
our boundary value problem can be found from a geometrical consideration. First of all, 
we note that the function 

, . , ln|0'-0| 
g{g>,0):=-L- 1 (A6) 

Z7T 
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satisfies AQ(0 l , 0) = 8(0 — 0'). Hence, the vector field 



Ju A 0-0'r Ju 



satisfies the 'Poisson' equation (Al), i.e., V • h(0' , 0) = -5(0 - 0') + 1/A, with A = L ± L 2 . 
The vector field h(0',0) has the geometric interpretation of a deflection angle at point 
0' caused by a homogeneous mass distribution of surface mass density 1/(2 A) inside the 
rectangle, and zero surface mass density outside, plus that of a point 'mass' at position 
in U with negative mass —1/2. Of course, h does not satisfy the boundary condition (A2). 

Let us define the points {ij) for i ^ 0, j ^ by {ij) = (flj , 0^) , with (k = 1, 2) 

0(0 = J + (* ~ l ) L k for i odd , ,^ g * 
k \ iLk — 9k for i even, 

for i > 1, and with = — 0^ for negative i. Furthermore, we define the rectangular 
regions to have sidelengths Li, L 2 , and to be centered on = (ci^> c 2^)> with 

4° = (2i - l)L k /2 for i>\ , 
cj^ = —c k for z < — 1 

Then, (11) = 0, W( n ) = U, 0^~^ = -0^\ etc. Note that the point {ij) lies in 
Ui % i) . Fig. 11 illustrates the geometrical arrangement of the points 0^ and the rectangles 
US % ^ . In analogy to (A7), we define the vector fields, for i ^ 0, j ^ 0, 

h^(0',0) = -vg(0 l ,0^A+ [ d 2 ^ vg ^ 0/ ) , (^9) 

which can be interpreted as the deflection angle at the point 0' caused by a uniform surface 
mass density of magnitude 1/(2 A) in the rectangle plus the deflection by a negative 
point mass -1/2 at {lj) . We then have, for e U, 0' e U, that h (11) (0', 0) = h(0' , 0), and 
V • h^(0', 0) = for (i, j) ^ (1, 1). Since the geometrical arrangement of rectangles and 
points is symmetric with respect to a reflection around the lines 0[ = mLi, and 6' 2 = nL 2 , 
the sum 

u(0', 0) = J2 E hW) (^) 

1*1=1 b\=i 

satisfies the Poisson equation V-H(0', 0) = —5(0 — 0') + l/A for 0, 0' e 14, and also satisfies 
the boundary condition (A2), due to symmetry, provided the sum in (A10) is extended to 
infinity in all directions. However, special care has to be taken by performing this sum: 
the magnitude of decreases as 1/r 2 , where r is the distance of the rectangle from 
the origin (and not as 1/r, since the 'monopole' contribution of the negative point mass 
is canceled by the uniform surface mass density, i.e., the total mass inside vanishes), 
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Fig. 11. The arrangement of the fields 'US 1 ^ and angles [ s scetched for —4 < i, j < 4, i,j ^ 0. The 

stars denote the positions c( IJ ) = i^c\^ , c ^^j of the centers of a fields U^ % ^\ the filled squares denote the 

points Q^ 1 ^. The boundary of the field lA^ 11 ^ = U is marked with a thick solid line, the dashed lines are 
the coordinate axes 

and the two-dimensional sum can diverge for N x ,N y — >■ oo if this limiting process is not 
defined carefully. 

In order to see how this limiting process should be done (in particular for the numerical 
realization), consider the contribution of the four rectangles around the origin, 

w(0',0):= h W) (0',0) = - V0(0,0 W) ) + F(0',Li,L2) , (AH) 

i,j=±i 4i=±i 

where 

r A r B x - x' 

F(x\A,B)= / dx[ dx' 2 ^ . (A\2a) 

J-A J-B |x — x'| 

Evaluating the integrals, this becomes 
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Fx(x; A, B) = B + X2 In + ^) 2 + (A + xQ 2 + g - x 2 ^ (B - x 2 ) 2 + (A + xi) 2 



+ (A + Xl ) 



(A- Xl ) 



(B + x 2 ) 2 + (A - xi) 2 
B + x 2 



arctan 



arctan 



A + Xl 

B + x 2 
A - xi 



+ arctan 



+ arctan 



2 ~ (b - X2 y + (a - Xl y 

( B-x 2 
\A + Xl 
B-x 2 



F 2 ( X -A,B) = ^ln[ A + X ^ + ^ + X2) l + A 



(A + Xl ) 2 + (B-x 2 y 



A - xi 

- Xl (A - Xl ) 2 + (B + x 2 ) 2 
m 



+ (B + x 2 ) 
~(B- x 2 ) 



arctan 



arctan 



4±^U arctan 
B + x 2 J 

A + x ± 

B-x 2 



+ arctan 



2 - (A - Xl ) 2 + (B - x 2 ) 2 
A — xi 
B + x 2 
A - xi 
B-x 2 



We note the following expansion for F: 

tw , ^ f B \ AAB 

Fi( X , A,B) = 4 arctan Xl + 3( ^ 2 + ^ 2)2 



(Al2b) 



xi (x\ - 3x 2 2 ) 



4AB(A 2 -B 2 ) , A 9 9 A s ,„ 



5(A 2 + B 2 ) 1 



A 



A 6 



F 2 (x, A,B) = A arctan ( — ) x 2 + 



4AB 



3(A 2 + B 2 )' 



;X 2 (xl - 3xj) 



4AB(A 2 - B 2 ) 
5(A 2 + B 2 ) 4 



x 2 (hx\ - \§x\x\ + x\) + O 



A 6 



AABxi 



1 + 



A 2 



(Al2c) 



3A 4 - 10A 2 B 2 + 3B 4 xi - lOxlx 2 . + hx\ , n 
+ 1 -V +° 



A e 



F 2 ( X ,A,B) 



AABx 2 
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1 + 



3A 4 - 10A 2 B 2 + 3B 4 hx\ - lOxlx 2 , + xl ,„ 
+ 1 — P + 



A 6 



15 



The leading order term of w(0' , 6) for \0'\ — >• 00 behaves like \0'\ , because the dipole of 
the mass distribution inside these four rectangles vanishes. We will term this collection of 
four rectangles (and four negative point masses) a 'big rectangle'. Consider now this big 
rectangle to be centered on a point c'*- 7 ' = (2iLi, 2jL 2 ); then its contribution to the sum 
in (A10) will be w (o' — O^J . Even better behaved for \i\ , \ j\ — >■ 00 is the sum of four 

such big rectangles, 
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defined for i, j > 1, since the quadrupole of these four big rectangles also vanishes. 
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Fig. 12. This figure shows the arrangement of the 'big rectangles': the long -dashed lines mark the axis of 
the coordinate system, the thick solid line the boundary of the field U = JT^ 1 ' 1 ) as in Fig. 11. The short- 
dashed lines incidate the 'small rectangles' the big rectangles are drawn with solid lines, where their 
centers are denoted by triangles on the a;2-axis, squares on the a;i-axis, and pentaga elsewhere. Four of 
the pentaga are filled, to indicate how the sum over 4 'big rectangles', defined in equation (A13), has to be 
taken. The values for i and j at the left and the top of the figures denote the indices of the 'big rectangles' 
and their centers C^-?) 

An efficient way to perform the sum in (A10) is thus to split it into an 'inner' and an 
'outer' part, 

H(0',0) = H in (0',0) + H out (0) , (AU) 

where 

K in (e',e) = -J2Y,^. ^ + ^- I F(0';N x L l ,N y L 2 ) , (Alb) 

H=i|j|=i |0 (y) - 0' 

and we have combined the contributions of the second term in (A9) from the inner part, 
and the function F was defined in (A12). We will choose N x and N y to be odd. Then, the 
outer part of H is evaluated as a sum over the W's, 
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H out (0',0) = - 



M v 



J2 w {io) (o',o)+ Yl w {oj) (o',o) 

i=(N„ + l)/2 j=(N y + l)/2 



+ 



M„ 



(N m -l)/2 M y Mx 

i=l j=(N y + l)/2 i=(N x +l)/2j = l 



(A16) 



Using (A12c), we can now expand the functions w'*- 7 ' in terms of 1/ 
2nW^\e', 6) = (q lVl + q 2 v 2 + q 3 v 3 ) ( "f 1 ) 
+ (10g 2 wi + 7q 3 v 2 ) 



, and obtain 



39?) 
6'P 



+ 21q 3 v 1 



-0[ (0? 

0' 2 {30? - 0?) 

-9[ {O'f - 1O0' 1 2 9' 2 2 + 50' 2 4 ) 
6' 2 (5^ 4 - 1O9 ,2 0' 2 + 9' 2 4 ) 



(A17) 



+ 



with 

V3 

Qi 

Q2 



L\ - 3(9 2 



el) , 

SiLi + L^-lO^L^ 2 -^(ei- 
SlLt-Ll-ri^^iLl-Lj)] 
16 (C 4 - 6C 2 C 2 + C 4 ) 
|C| 8 

16[Cj» - 15(dC 2 ) 2 (C 2 - C 2 ) 



60 2 2 ) 



2l[^-e 2 6 -15(^e 2 ) 2 ^ 2 -^ 2 )] 



(A18) 



CI] 



3|C 



12 



93 



16 [C? + C| - 28(dC 2 ) 2 (Cf + C 4 ) + 70(dC 2 y 



3|C 



16 



where we have written C = Since we are interested in the field H only for 6 and 

0' within the rectangle U, of size L = 0(Li, L 2 ), the next order term in (A17) is of order 
L 9 1 |C| 10 . We then note that only the coefficients q n depend on the center position 
hence, the sums of (A16) only need to be performed over the q n , whereas the v n depend only 
on 0. It is clear that this summation can be carried out very quickly, even if the number 
of terms (~ M x M y ) is large. We have evaluated the kernel H(0',0), using the preceding 
equations. Due to the explicit symmetry of our evaluation, the boundary condition (A2) is 
automatically satisfied at 0[ = and at 0' 2 = 0; at the other two sides of the rectangle, the 
boundary condition is satisfied only if a sufficient number of terms are taken into account. 
The value of H • n at the two remaining sides thus provides a measure for the accuracy of 
the calculation. It turns out that using N x = N y = 7 in (A15) and M x = M y = 400 in 
(A16) yields an accuracy of better than 10~ 6 . 

The kernel H was calculated on a regular grid, with 



0<i< Nx 



0<j<N 2 
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and 



w = ((fc-V2)£iM, -1/2)^2/^2) , l<k<N! , 1<1<N 2 . 

Of course, for a square it is useful to choose Ni = N 2 ; otherwise, the gridcells should be 
chosen to be nearly quadratic. The numerical routine is very efficient and has on output 
the two arrays Hi(0'-,Oki) and H 2 {0'-, Oki), i.e., arrays of dimension (JVi + 1) x (N 2 + 
1) x JVi x N 2 . For example the computation time for Ni = N 2 = 50 is about 20 minutes 
on an IBM rise 6000 processor. 
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Fig. 1. Contour- and surface plot for the function — _fC(x) = — In [1 — k(x)], where k is the two-dimensional 
surface mass density of the lens A,B,C and D respectively. In each panel, the size of the field is 7. '5 X 7/5, 
and K is calculated on a 50 X 50 grid; the spacing in the contour lines is 0.1 for lens A and lens B and 0.05 
for lens C and D 
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Fig. 10. Root mean square error of the reconstruction as a function of position for the KS-, SS-, S- and 
B-reconstruction. The number of simulations is 500, no = 50, p = 0.3, 6q = 0.3 and sm = 0.11; the 
corresponding power spectra for the error field are displayed in Fig. 5b. The vertical axis is drawn from 
to .2; the contour lines are in the same interval, with a step size of 0.0125; the size of the CCD field is 7.5 
arcminutes, as usual 
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